Dynamical Screening and Ferroelectric-type Excitonic Instability in Bilayer Graphene 
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Electron interactions in undoped bilayer graphene lead to instability of the gapless state, 'which- 
layer' symmetry breaking, and energy gap opening at the Dirac point. In contrast to single layer 
graphene, the bilayer system exhibits instability even for arbitrarily weak interaction. A controlled 
theory of this instability for realistic dynamically screened Coulomb interactions is developed, with 
full acount of dynamically generated ultraviolet cutoff. This leads to an energy gap that scales as a 
power law of the interaction strength, making the excitonic instability readily observable. 

PACS numbers: 



Graphene, due to its unique electronic structure of a 
two-dimensional semimetal, provides an entirelynew set- 
ting for investigating many-body phenomena Since 
the conduction and valence band joined together in a 
semimetal mimic massless Dirac particles, electronic phe- 
nomena in graphene often have direct analogs in high en- 
ergy physics [2|. In particular, several authors discussed 
the analogy between excitonic instability in a single-layer 
graphene and chiral symmetry breaking in 2-1-1 Quantum 
Electrodynamics. While this instability is absent when 
interactions are weakQ, and the situation for realistic in- 
teraction strength is still debated, the instability can be 
"catalyzed" by a magnetic field[l,[l|. These predictions 
are in qualitative agreement with experiment [6|. 

The effect of interactions is drastically different for 
semimetals with linear (type I) and quadratic (type II) 
electron dispersion. This was recognised in an early work 
0, where the difference in behavior was traced to the 
density of states at low energies, which is much lower in 
type I systems than in type II systems. Electronic prop- 
erties of quadratically dispersing systems are governed 
by infrared divergences in Feynman diagrams, resulting 
in unconventional low energy states, whereas in linearly 
dispersing systems the free-particle description is robust. 
Indeed, excitonic instability in BLG, which has gapless 
quadratically dispersing electronic states, was shown to 
occur for arbitrarily weak short range repulsion [8j] . 

In this Letter, we analyze excitonic instability in a 
BLG system with 1/r interaction, focusing on a fer- 
roelectric (FE) state that spontaneously breaks which- 
layer symmetry and polarizes the layers in charge. After 
accounting for dynamically generated ultraviolet cutoff, 
treated in the RPA screening framework, we find a gap 
which in the weak coupling limit scales as a square of the 
interaction strength, A c>c (c^/k)^, with e the electron 
charge and k, the dielectric constant. This is in contrast 
to the exponential BCS-like behavior of the gap expected 
in single layers coupled via a dielectric spacer [9- 11]. 

Interestingly, we find the behavior to be highly sensi- 
tive to the specifics of screening model: log divergent 
diagrams can become log^ divergent upon going from 
static to dynamic screening [see FiglUb)]. Thus for a 
reliable estimate of the gap it is necessary to properly 
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FIG. 1: (a) Domains of opposite polarization in the ferroelec- 
tric state. Valley polarized chiral edge states propagate in 
opposite directions along domain boundaries, (b) Diagram- 
matic representation of gap equation. First term is vertex 
correction, second term is self energy correction. Both dia- 
grams exhibit log^ divergence which cancel to leading order. 
Solid, wavy and thick blue lines represent fermion propaga- 
tors, the RPA interaction ([4]), and Ars vertex, respectively. 

treat dynamic screening. We evaluate the dynamical po- 
larization function for BLG and use it to estimate the 
gap value. For realistic parameters the predicted gap lies 
within the experimentally accessible range. Generaliz- 
ing our approach for the anti-ferromagnetic (AF) states 
considered in we find that the AF and FE states are 
nearly degenerate in energy. 

The formation of a gapped state will manifest itself in 
strongly temperature dependent conductivity at T < A. 
In the presence of long-range disorder, the gapped state 
will occur at the p-n boundaries separating electron and 
hole 'puddles', making these boundaries a bottleneck for 
transport. Hopping- like temperature dependence is in- 
deed noted in all conductivity measurements near neu- 
trality point in BLG but not in single-layer graphene [l[. 

Since the ferroelectric transition breaks a Z2 symme- 
try, excitonic ordering will produce domains of opposite 
polarization [see FiglUa)]. In the absence of disorder, 
the characteristic size of the domains L is determined 
by long-range attraction between polarization in neigh- 
boring domains, in analogy with "electronic microemul- 
sion phases" discussed in Ref. 17|, giving an estimate 
8p^ln(L/ao) ~ A/ao, where p is polarization density 
and ao is the correlation length, Eq.®. As discussed 
m Ref.fl^, the boundaries between regions with oppo- 
site polarization host valley polarized edge states. Since 
a two-dimensional system with two domain types should 
exhibit percolation of edge states, the FE state should 
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be able to carry valley currents. This could be useful for 
valley tronics applications. 

The low-energy Hamiltonian for BLG can be described 
in a 'two band' approximation, neglecting the higher 
bands that are separated from the Dirac point by an 
energy gap W ~ 0.4 eV [12|. The electron states are 



described by wavefunction taking values on the A and 
B sublattice of the upper and lower layer respectively. 
The non-interacting spectrum consists of quadratically 
dispersing quasiparticle bands E± = ±p^/2m with band 
mass m ~ 0.054mp. It is convenient to introduce the 



Pauli matrices that act on the sublattice space, and to 

Px±iCPy UM' where C 



space 

define t± = ti ± iT2 and p± = px i: iC,Py 
= 1 for the K valley and ^ = — 1 for the K' valley. The 
Hamiltonian may then be written as 



H 



Ev4 



2m 



(1) 



Hn 



n(x)n(x') 



2k ^ Ix 



The sum over a indicates summation over TV = 4 spin 
and valley species, while the dielectric constant k incor- 
porates the effect of polarization of the substrate and 
of the higher bands of BLG. The interaction is invariant 
under SU{N) rotations in spin/valley space. We also ap- 
proximate by treating the interlayer and intra-layer in- 
teraction as equal, and defer discussing the effect of finite 
layer separation until after we present our main result. 

We investigate stability of the gapless state by intro- 
ducing a test gap-opening perturbation At^ into the non- 
interacting Hamiltonian, where A must be real, but may 
take either sign. This test perturbation explicitly breaks 
the Z2 layer symmetry of the Hamiltonian, and corre- 
sponds to a ferroelectric instability that polarizes the lay- 
ers by charge. We develop our analysis perturbatively in 
the interaction, and calculate the interaction renormal- 
ization of the Ara vertex. At leading order in weak bare 
interactions, the vertex correction in Figlljb) takes the 
form ST — tsSA, where 



(lE + Ho){iE - Ho 



-U. 



(2) 



Here Hq is the Hamiltonian of the non-interacting sys- 
tem evaluated at A = 0. The vertex correction is posi- 
tive and preserves the form of the (73 vertex. Moreover, 
simple power counting shows that the vertex correction 
is divergent in the infrared for any form of interaction 
U, screened or unscreened. The infrared divergence in- 
dicates instability even for arbitrarily weak interactions 
(unlike monolayer graphene). The infrared divergence is 
power law when U is the unscreened Coulomb interac- 
tion, therefore it is important to include screening even 
at weak coupling, to moderate the infrared divergence. 



We now introduce the interaction energy scale Eq and 
the corresponding lengthscale qq, defined as 



Eo 



1.47 



eV, 



ao 



K X 1.1 nm. (3) 



For simplicity, we take the weak coupling limit Eq < W, 
and neglect interaction-induced mixing of the low energy 
states with the higher bands [for discussion of mixing 
see Ref.fl3|]. Moreover, when U is the unscreened or 
dynamically screened Coulomb interaction, the integral 
in Eq.(I2]) is convergent in the ultraviolet limit, without 
the need for any high energy cutoff. Thus Eq emerges as 
the only energy scale in the problem. This then implies 
that the energy scale for the gap must scale as a power 
law in electric charge A ~ i?o ~ e^. 

As we shall see, it is necessary to properly treat dy- 
namic screening to obtain a reliable estimate for the gap. 
We therefore take U in Eq.Q to be the dynamically 
screened Coulomb interaction, defined as 



2ne' 



Kq - 



2Tre^NU, 



(4) 



a;,q 



Here we have introduced the single species polarization 

, ,2 

function Htj,q = — / G{p+,e+)G{p^,e^)-^^, where we 
use the notation p± = p ± q/2 and e± = e ztuj/2, and 
define the imaginary frequency Green function in terms of 
the non-interacting Hamiltonian Hq as G^^{E,p, A) = 
iE — Hq{p, A). Hence, we obtain 



H, 



= -2 



dE(Pp E+E- - pIp± cos(26'pq) - A 
(2^)3 {El + pX + A^){El 



-P- 



A2)- 

(5) 

Here, 9pq is the angle between the vectors p+ and p_. 
To determine the dynamically screened interaction, it is 
sufficient to determine the polarization function in the 
ungapped state. We therefore suppress A in Eq.Q, in- 
tegrate over frequencies by residues, and scale out q. The 
integral then depends on a single dimensionless param- 
eter oj = 2mu!/q'^, and may be evaluated analytically 
by integrating over momenta in polar co-ordinates (see 
Appendix) . This gives an exact expression for the polar- 
ization function H^j^q^o = ~|^/('^)j where 



tan 2u} 



In 



a;2 



1 
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(6) 

where u = (21n47r)2. The right hand side provides an 
approximate formula that reproduces /(w) exactly for 
cj — > and w — )■ 00, interpolating accurately between the 
two limits. Result ([6]) agrees with the polarization func- 
tion found in [isj continued to Matsubara frequencies. 

The vertex correction, Eq.®, calculated with the dy- 
namically screened interaction, Eq.Q, is log^ divergent 
in the infrared (see Appendix). The enhanced divergence 
arises from the phase space region uj/q'^ ^ 1, where the 
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a) 





c) 



V 



d) 



FIG. 2: Free energy change from gap formation at leading or- 
der in A and in the interaction (notation the same as in Fig[T} . 
The diagrams may be interpreted as: a) Kinetic energy cost 
from spontaneous gap opening; b, c) Interaction energy gain 
from gap opening; d) Hartree energy cost of layer polarization, 
which vanishes within the approximations of Eq.([T]). Here 
V- is the difference between interlayer and intra-layer inter- 
actions. While all these diagrams are nominally O(A^), A 
also appears as a logarithmic infrared cutoff in each diagram, 
Eq.®. 



the RPA approach, by 

/I 72 
-^-^ In (l - iVi7„,g(n^,,,A - n^,g,o)) . (7) 

Here, U^.q is the interaction (U) and II^.g.A is the single 
species polarization function in the gapped state. The 
problem thus reduces to that of evaluating the polariza- 
tion function at finite A. 

We calculate the quantity Ha — Hq by integrating 
Eq.(l5]) over frequencies by residues, Taylor expanding to 
leading order in small A^, and analytically performing 
the integration over momenta, assuming as before that 
A is independent of momentum. After some algebra (see 
Appendix), we obtain 



27rr2 



-7 



,9 



In- 



(8) 



Coulomb interaction is not efficiently screened. However, 
it is also necessary to take account of the self-energy cor- 
rection, so that the full gap equation is represented by 
FiglUb). In particular, the self energy undergoes a log^ 
renormalization 2l| , and this can be shown to cancel the 
vertex correction at log^ order, leaving a residual loga- 
rithmic divergence. 

Since demonstrating the cancellation of the self energy 
and vertex correction at leading log^ order and extracting 
the subleading logarithmic divergence is fairly tedious, we 
employ an alternative scheme for solving the gap equa- 
tion Figlljb). We note that calculating the free energy 
of BLG as a function of A, at leading non- vanishing or- 
der in A, simply produces the diagrams in Figl2lja,b,c). 
Upon minimizing with respect to A, this yields the gap 
equation, pictured in Fig [If b), with correct combinato- 
rial coefhcients. Minimizing the free energy of BLG with 
respect to A is therefore formally equivalent to solving 
the gap equation, and is technically simpler. It may be 
verified that while FigH] (b) and (c) are individually log"^ 
divergent in the infrared, their sum is only log^ divergent. 
This is the same leading order cancellation of divergences 
that is manifested by the gap equation, Figlljb). 

We approximate by assuming that the gap function A 
is static and momentum independent up to energies of 
order Sq, on the grounds that the screened interaction 
in the particle-hole channel depends only weakly on the 
transferred momentum. We evaluate the kinetic energy 
change 5T represented by Figl2ja) by including A in the 
fermion Green function. We find, with logarithmic accu- 
racy, 6T= ^A"^ ln(A/A), with the cutoff A^^ Eq. 

To calculate the exchange energy gain, we note that 
the difference in interaction energy between the gapped 
and ungapped states, SE = E{A) — E{0) is given, within 



We are using the notation — uj^ + q'^. 

We now evaluate the exchange energy gain from gap 
formation by substituting Eq.® into Eq.([7]), and per- 
forming the integrals using polar coordinates (see Ap- 
pendix). Combining this with the kinetic energy cost 
ST, we find the free energy associated with gap opening 



^^(A)==-A2ln(A/A) 



13toA2 



\ii^{N^Eq/A). (9) 



We note the expected emergence of a natural ultraviolet 
cutoff. Identifying A with N'^Eq, and minimizing Eq.([9]) 
with respect to A, we obtain, with logarithmic accuracy. 



A = N^EQexp{-3TT^N/13). 



(10) 



We emphasize that Eq appears only outside the expo- 
nential, making A a power law function of interaction 
strength at weak coupling. However, A is exponen- 
tially small in N, where N is the number of fermion 
species participating in screening. If we had worked 
instead with static screening, we would have obtained 
A ~ iV^E'o exp(— 2A^ln4), and would have underesti- 
mated the size of the gap by an order of magnitude. 

For iV = 4 Eq.ldni) gives A ^ IO'^Eq « LS^-^meV, 
upto a numerical prefactor of order unity. Meanwhile, 
numerically evaluating the integrals in Eq.([71) and the 
kinetic energy contribution ST [Figl2][a)], and minimizing 
the free energy gives A « 4K~'^meV. This number lies 
within experimentally accessible range, although it can 
be reduced by screening in the substrate, by doping, or 
by disorder induced density inhomogeneity. 

Thus far, we have neglected the effect of trigonal warp- 
ing which leads to deviation of particle dispersion from a 
simple quadratic dependence, causing an overlap of the 
conduction and valence bands. Trigonal warping can pro- 
vide an alternative low energy cutoff, preventing forma- 
tion of a gapped state if its size exceeds the estimated 
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gap. Since the the upper estimate for trigonal warp- 
ing, 1.5 meV [2^, is less than the numerical estimate 
A « 4K~^meV, the effect of trigonal warping should be 
inessential, at least for suspended bilayers (k « 1). 

Our analysis can be easily generalized to any state that 
adds a term Arsfi to the Hamiltonian ([T|), where 17 is a 
4x4 hermitian matrix in spin/valley space satisfying 
D,'^ — 1. The FE state considered above corresponds to 

— 1, wheras the AF states discussed in Ref.Q are char- 
acterized by = 0-3(8)1 or = l(g)ri3, where as and 773 are 
Pauli matrices in spin and valley space, respectively. All 
these inequivalent choices for fl yield the same mean field 
free energy F{A), Eq.®, and the same gap value as was 
obtained for the FE state. This mean field degeneracy 
occurs because the Hamiltonian is invariant under SU (4) 
spin/valley rotations, within validity of Eq.([T]), while the 
states corresponding to different choices of il differ only 
in their spin/ valley structure. 

We now examine the effect of finite layer separation 
d ~ SA, which differentiates the interlayer and intra- 
layer interactions, giving an anisotropy V- — ^{Vaa ~ 



Vab) 



1 27re 

2 q 



^(1 - e-""*) = 



ire'^d. This anisotropy is 



small, because d <^ qq. The leading order effect of finite 
layer separation is to introduce a Hartree energy for the 
states that polarize the layers in charge [see Figl2][d)], 

Afm^ 

i^Hartrcc = -T^V- \n\A/A). (11) 

It was argued in Ref.poj that this contribution prevents 
formation of the ferroelectric state. However, Ref.[20l| 
neglected the exchange energy. We note that the Hartree 
energy is of the same functional form as the exchange 
energy, Eq.®, but is parametrically smaller by d/ao ^ 
1, and so cannot prevent the instability. 

Upon going beyond the weak coupling approximation, 
d/ao ceases to be a good control parameter, but our con- 
clusions remain unchanged. This is because the in- 
teraction is screened as V- — V^/{1 ~ V^_H_), where 

H_ = /TrG(e+,p+)T3G(£_,p_)r3f$ ~ ^WA). 
Such logarithmic screening ensures that the Hartree en- 
ergy remains smaller than the exchange energy, Eq.®, 
and so cannot prevent gap formation. However, when 
d/ao is not small, the Hartree energy may tip the bal- 
ance from the FE state to one of the AF states, which do 
not polarize the layers by charge. 

A ferromagnetic instability was predicted for un- 
screened interactions in [l^]. However, the free energy 
gain from ferromagnetism was only cubic in the ferromag- 
netic order parameter. In contrast, the excitonic states 
all have a free energy gain of O(A^), and should thus 
dominate for weak coupling. 

After submitting this Letter we became aware of works 
22 [ and Instabilities in BLG are analyzed in these 
papers within a renormalization group framework, with 
interactions modeled as being short ranged. It is found 



that different choices of short range interaction can result 
in different states, gapped [2^ or gapless 23 1. 

To conclude, we have demonstrated that electron in- 
teractions drive BLG to a gapped state. For dynamically 
screened Coulomb interactions, taking proper account of 
dynamically generated ultraviolet cutoff and log^ diver- 
gences arising in this case, we obtain a gap value in ex- 
perimentally accessible range (A « 4K~^meV). Mani- 
festations of the energy gap opening in BLG at charge 
neutrality include temperature dependent transport at 
low temperatures, and also valley polarized chiral edge 
states propagating along domain boundaries. 

We acknowledge useful discussions with D. Abanin, 
P. A. Lee, Nan Gu, D. Mross, and A. Potter. This work 
was supported by Office of Naval Research Grant No. 
N00014-09-1-0724. 
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APPENDIX A: CALCULATING POLARIZATION 
FUNCTION 

Throughout this supplement, we use a shorthand 
whereby represents /2m, unless otherwise specified. 
We wish to evaluate the polarisation function at non-zero 
A. This was defined as 



n, 



aj,q,A 



dEcPp E+E^ - p^pt cos(2gp^) - A 



(27r)3 {El +p\ + A2)(£;2 +pt+ A2) ' 

(12) 

where E± = E ±oj/2, p± = p ± q/2 and 6pq is the angle 
between p+ and p_ . Wc begin by integrating by residues 
over frequencies, to obtain 



n, 



= - [ ^ 

J (27r: 



)2w2. 



2sin{9pgf 



~2 ~2 

+ cos(2^p,)(l-^) 



We are using the notation ~ p^. + A"^. Upon set- 
ting A = 0, we obtain the polarisation function in the 
ungapped state. We take q to lie along the x axis 
without loss of generality, change to polar coordinates 
Px = p cos tp, Py = p sin If and scale out q. The inte- 
gral then depends on a single dimensionless parameter 
uj = oj/q'^, and may be evaluated analytically by inte- 
grating in turn over ip and p/q. This gives an exact ex- 
pression for the polarisation function n^j.q.o = —^f{^)^ 
where 



2 tan" 



CO — tan 



^2(1; 



+ ln- 



^2 + 1 



In 4 



(14) 

with u= {2 ln4/7r)2. We note that the polarisation fimc- 
tion vanishes for g2 <^ go that in this regime of pa- 
rameter space, the interaction is not efficiently screened. 



We now calculate IIa — Ho by Taylor expanding IIa to 
(13) leading order in small A2 to obtain 

_J 



Ha-Ho = A' 



d'^p 



pI+P- 



f 1 



(27r)2a;2 + (p^+p2_)2 ^^^^2 



\{\ + \)cOS2epq 

2 p% pt 



sm 



+ - 



p+p't , 



(u;2 



{Pl 



.p2_)2)2 



(15) 



We note that the integral has singularities at p± = 0, 
which must be regularised by introducing an IR cutoff 
A. The singularities at p+ and p_ contribute equally, so 
we evaluate just the contribution from the singularities 
at p- = 0, and multiply by two. We change variables to 
p_ p, and p+ — >^ p + q and work in polar coordinates, 
taking q to lie along the x axis. Again, we scale out q 
and keep only the singular terms, to obtain 



power series in small p. We note that the terms at 0(p~^) 
and 0{p~^) vanish upon angular integration, so that the 
leading term comes at 0(p~2) After integration over 
angles if, and after making the substitution z = p^ /2m 
we obtain 



Ota = 



mA2 



-7 



+ 



1 



ATT ^X2 [ Pdpdy 1 1 , , 
dllA = 2A^ / , ,„ — 5 — -^H , „ — -F{p) 

J (27r)2w2 + lp2 (w2 + l)2 2 ^' 

(16) 

Here, we have measured all variables in units of q, have 
defined ip to be the angle between p and q, and have 
defined 



2-K \{uJ^^ + lf (a;2 + l4)3 
After restoring q, we obtain 



/fas, 



mA2 iV^±7. 



+ 



4g8 



(a;2 + qi)2 {^2 + ^4)3 ) ■ 
(19) 



F{p) 



1 + 2p2 + 2p cos (/? 1-2 sin2 6pq 

up- -h (2j32 + 2pC0S</3-|- 1)2 ^4 



(17) 



We can readily determine that in this basis, sin B^q = 
sin2 (p/(p2 -(-2pcos(/5-|- 1). We can then expand -F(p) as a 



APPENDIX B: CALCULATING EXCHANGE 
ENERGY 

In the main text, we derived 

SE = J -^-^ In (1 - 7vc/„,,(n„,,,A - n„,g,o)) . (20) 
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Substituting the result p9|) into this, and expanding the logarithm to leading order in small A^, we obtain 



exchange 



(27r)3 27r 



ln(- 



A 



4q8 



27r 



I 



(21) 



We have used the approximate form for the interaction, 
and have measured everything in units of Eq and oq. We 
note that Ha — Hq vanishes in the limit /uj <^ 1, when 
the Coulomb interaction is not efficiently screened. This 
is the region of phase space that is responsible for the 
vertex correction having an enhanced log^ infrared di- 
vergence (see next section). The vanishing of IIa — Ho h^ 
the region of phase space where the Coulomb interaction 
is strongest is responsible for the weakened divergence, 
and is another way of understanding the leading order 
cancellation between the vertex correction and the self 
energy. To formally calculate the exchange energy, we 
approximate the interaction by taking 



1 2d{N - q) 

q + Nq^hiA/ ^q-^ + {2\nA/ nfuj'^ Niiq^/uj 



(22) 



where 9{x) represents the Heaviside step function. We 
take the q^ > to limit of the interaction, because this is 
where the interaction is strongest. We now substitute 
into Eqn.lO, change to polar co-ordinates uj = rcosO, 
q^ = rsin^, and integrate in turn over 9 and r to obtain 
the result quoted in the main text 



E, 



exchange 



-^\n\N^Eo/^) (23) 



APPENDIX C: THE log^ DIVERGENCE OF 
VERTEX CORRECTION 

We analyse the vertex correction 5A with the dynami- 
cally screened effective interaction. After including A in 
the fermion green functions, we obtain 



dedzdO £ 



A2 



27re2A, 



(2^)3 (£2 



z2 ^ A2 p)2 k|p ■ 



me'^Nf{uj) 
(24) 



We use the notation z = and ui — uj / \p — q\'^ . We 
now measure energies and wavevectors in the units Eq 
and l/flQ respectively, to recast the vertex correction in 
the dimensionless form 



dedzdO 



£2 + z2 



A2 

e.p 



27rA^ 
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(25) 

where z' = (p — q)2 and ^ = (£ — '-^)/z' . We note that 
integration over z can be extended to infinity without 
running into a UV divergence. 

We approximate by treating the gap function A as 
a constant. Using the approximate form of /(w) ~ 
ln4/vlT?«2?, with u = (21n4/7r)2, and changing to 
polar coordinates e = p&inip^ z = pcostp, we integrate 
in turn over p and cp. Working to leading order in small 
A/{N'^Eo), we obtain 
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One logarithm comes from the integration over 'an- 
gles' (fi, and rcfiects the singularity of the dynamically 
screened interaction in the region q^ /u ^ 1, whereas the 
second logarithm comes form the integration over p, and 
reflects the non- vanishing density of states at low energies 
in BLG. This second logarithm would be present even if 
we worked with short range (delta function) interations. 

The self energy is also log^ divergent, for the same rea- 
sons, and cancels the logarithmic divergence that comes 
from the singularity of the interaction as g — ?► 0. This 
manifests itself in the free energy calculation through 
IIa — Ho vanishing as ~ in the limit q 0. However, 
we note that the cancellation between vertex correction 
and self energy exists only to leading order. At sublead- 
ing, logarithmic order, there remains an instability, as 
was made evident through the free energy approach. 



